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1. Introduction 

Over the last few years the detection of several giant air showers, either through ground 
based detectors [TJ |2] or fluorescence telescopes [3 E] , have confirmed the arrival of 
ultra high energy cosmic-rays (UHECRs) with energies up to a few hundred EeV (1 
EeV = 10 18 eV). Their existence poses a serious challenge and is currently subject of 
much theoretical research as well as experimental efforts (for recent reviews see jSlElIZI)- 

The problems encountered in trying to explain UHECRs in terms of "bottom- 
up" acceleration mechanisms have been well-documented in a number of studies (e.g., 
Refs. [HJ El EH!)- I n summary, apart from energy draining interactions in the source, 
the maximal UHECR energy is limited by the product of the accelerator size and the 
strength of the magnetic field. According to this criterion it turns out that it is very 
hard to accelerate protons and heavy nuclei up to the observed energies, even for the 
most powerful astrophysical objects such as radio galaxies and active galactic nuclei. 

In addition, nucleons above ~ 50 EeV suffer heavy energy losses due to photo- 
pion production on the cosmic microwave background (CMB) — the Greisen-Zatsepin- 
Kuzmin (GZK) effect jib - which limits the distance to possible sources to less than 
~ 100 Mpc [12J. Heavy nuclei at these energies are photo-disintegrated in the CMB 
within a few Mpc [T3]. Unless the sources are strongly clustered in our local cosmic 
environment, a drop, often called the "GZK cut-off" in the spectrum above ~ 50 EeV 
is therefore expected [Tl] . even if injection spectra extend to much higher energies. 
However, the existence of the latter is not established yet from the observations |15j . In 
fact, whereas a cut-off seems consistent with the few events above 10 20 eV recorded by 
the fluorescence detector HiRes jl], it is not compatible with the 8 events (also above 
10 20 eV) measured by the AGASA ground array [2]. The solution of this problem may 
have to await the completion of the Pierre Auger project JH] which will combine the 
two complementary detection techniques adopted by the aforementioned experiments. 

Adding to the problem, there are no obvious astronomical counterparts to the 
detected UHECR events within ~ 100 Mpc of the Earth HE Ej. At the same time, 
no significant large-scale anisotropy has been observed in UHECR arrival directions 
above ~ 10 18 eV, whereas there are strong hints for small-scale clustering: The AGASA 
experiment has observed four doublets and one triplet within 2.5° out of a total of 
57 events detected above 40 EeV [2|. When combined with three other ground array 
experiments, these numbers increase to at least eight doublets and two triplets within 
4° [15] . This clustering has a chance probability of less than 1% in the case of an 
isotropic distribution. 

There are currently two possible explanations of these experimental findings. In the 
first one, which assumes negligible magnetic deflection, most of the sources would have to 
be at cosmological distances which would explain the absence of nearby counterparts and 
the apparent isotropy would indicate that many sources contribute to the observed flux, 
where a subset of specially powerful sources would explain the small-scale clustering [19]. 
This scenario predicts the confirmation of a GZK cutoff. The second scenario is more 
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realistic and takes into account the likely existence of large scale intervening magnetic 
fields correlated with the large scale galaxy distribution. In this case magnetic deflection 
could be considerable even at the highest energies and the observed UHECR flux 
could be dominated by relatively few sources within about 100 Mpc. Here, large scale 
isotropy could be explained by considerable angular deflection leading to diffusion up 
to almost the highest energies and magnetic lensing [20 could contribute to the small 
scale clustering. 

A possible way of explaining the origin of these particles is to assume a 
source population associated with quasar remnants currently inactive in the visible 
spectrum [21]. If these sources are sufficiently numerous, this would explain large 
scale isotropy, whereas a few nearby sources could explain the small scale clustering. 
In this scenario, UHECRs would be accelerated at the nuclei of those nearby dead 
quasars, whose underlying super-massive black holes are sufficiently spun-up to provide 
the necessary electromotive force. Hints of a possible correlation between some of these 
dormant quasars and the arrival direction of cosmic rays above 4 x 10 19 eV may support 
this scenario [22 . 

In the present study we perform numerical simulations in order to test a sample 
of candidate dead quasars as possible sources for the ultra high energy events. We 
first considered the case of an all-pervading turbulent magnetic field of InG strength, 
which corresponds to the primordial field necessary to explain galactic magnetic fields 
by adiabatic compression [22] ■ We also performed simulations with a relatively strong 
magnetic field of 0.1/iG, either following a pancake like structure mimicking the Local 
Supercluster, or homogenous up to the farthest sources. All these cases are consistent 
with Faraday rotation upper limits [21], although the strong field case may not be very 
realistic. 

The paper is organised as follows: in Sec. II we outline the scenario underlying our 
numerical simulations which are briefly described in Sec. III. In Sec. IV we introduce the 
statistical quantities used for comparison with the data. In Sec.V we present our results 
and we conclude in Sec. VI. 

2. The Quasar Remnant Scenario 

Boldt and Ghosh [21] suggested that supermassive black holes, present in the center of 
normal galaxies [22], could be the principal sources of particles at the highest energies. 
As stressed in the introduction, the lack of observed counterparts to the highest energy 
events within an acceptable distance suggests the idea of sources currently inactive in 
the visible spectrum. These objects should be able to accelerate UHECRs and appear 
dormant at the same time. Quasars remnants or dormant AGNs (active in earlier 
phases), with underlying supermassive black holes, are among the astronomical objects 
which best satisfy these requirements. In this model the black holes were spun up to 
nearly their maximal spins during an earlier phase, when the dormant AGNs were active, 
and liberate their rotational energy in the form of UHE particles rather than in powerful 
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jets. 

The acceleration is due to the potential difference produced by magnetic field 
lines threading the event horizon of the supermassive black hole in rotation. The 
electromotive force generated by a rotating black hole of mass M, threaded by a 
magnetic field of strength B extending on a range R ~ R g (gravitational radius 
= GM/c 2 ) is given by [2E] 

emf « 9 x 10 20 (a/M) B A M 9 volts (1) 

where a is the black hole specific angular momentum, B^ = B /(10 4 G) and M$ = 
M/(10 9 M Q ). For a rapidly spinning hole, we can approximate a/M ~ 1 which gives the 
maximal estimated energy. 

Apart from the mass of the object, the only input to estimate the maximal energy 
achievable in a given supermassive black hole is the magnetic field which will depend on 
the mass accretion rate. By assuming pressure equilibrium between the magnetic field 
and the infalling matter we have [21] 

B 4 = 1.33M^M 1/2 , (2) 

where M is the mass accretion rate in units of M Q /yr. From Eqs. Hand El we obtain 

emf « 1.2 x 10 21 M 1/2 volts. (3) 

In the present case, we consider protons as primary for which the energy loss dominant 
processes during the acceleration phase are due to pair production and photomeson 
production on ambient photons. We will neglect these energy losses at the source because 
an exceptional low accretion luminosity characterizes the supermassive black holes at 
the centers of nearby bright galaxies [21] making the mean free path for this interaction 
larger than R g . The losses due to curvature radiation induced by the magnetic field 
are taken into acount as in Ref. [2E]- Furthermore, as in [221, we assume the accretion 
rate onto the black hole to be about 10% of the bulge mass loss rate, M ~ O.IM^, 
with M± 2 = Mb u ig e / (10 12 M Q ) , and Mb u i ge the bulge mass. The maximal proton energy 
expected is then given by 

£ max = 137(M) 1/8 (M 9 ) 1 / 4 EeV 

~ 104(Mf 2 ) 1 / 8 (M 9 ) 1 / 4 EeV. (4) 

The list of sources we will use in our numerical simulations are given in Table Q 
and is obtained by selecting from the table of massive dark objects in nearby galaxies 
given by Magorrian et al. |25[ I2*H] . Kormendy (SOj, and Marconi et al. [2] those which 
could accelerate particles beyond 4 x 10 19 eV. We do not include NGC 4486 as from 
recent estimates of the spectral energy distribution at the core of this galaxy it is not 
possible to accelerate particles to the highest energy by the compact dynamo considered 
here [21] (second reference). Injection rates of individual sources are represented by a 
weight factor given by the mass accretion rate in units of M & /yr divided by the maximal 
energy in units of EeV. These weight factors are shown in the last column of Table ^ 
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Table 1. The list of sources used in our simulations. Here, D is the distance to the 
source in Mpc, M 9 = M/(10 9 M Q ), Mf 2 = M bulgo /(10 12 M Q ), E™ ax = E max / (10 19 eV), 
a and 5 are the equatorial coordinates and WF is the weight factor xlO 4 . 
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3. Numerical simulations 

The list of sources given in Table Q has been implemented into the numerical code 
for UHECR propagation in extragalactic magnetic fields used in earlier studies, see 
Refs. jH2l EH] for details. 

For each configuration many nucleon trajectories originating from the sources were 
computed numerically by solving the equation of motion for the Lorentz force and 
checking for pion production every fraction of a Mpc according to the total interaction 
rate with the CMB. In case of an interaction, secondary energies were randomly selected 
according to the differential cross section. Pair production by protons is treated as a 
continuous energy loss process. 

Each trajectory is abandoned if the particle reaches a distance from the observer 
twice bigger than the distance to the farthest source, i.e. 230 Mpc, or if the propagation 
time exceeds 10 Gyr. 

A detection event was registered and its arrival direction and energy recorded each 
time the trajectory of the propagating particle crossed a sphere of radius 1 Mpc around 
the observer. For each configuration this was done until 5000 events where registered. 
For each detected event, we register also the source by which it has been emitted. 

We assume a mono-energetic injection spectrum at the source: all particles coming 
from a source are emitted with the maximal energy of acceleration for that source as 
given in Table H We take into account the individual source power by including the 
weight factor in the last column of Table H 

We assume a random turbulent magnetic field with power spectrum < B(k) 2 >oc 
k nB for (2vr/lMpc) < k < (2?r/0.01Mpc) and (B(k) 2 ) = otherwise. The magnetic field 
modes are computed on a linear grid in momentum space and are Fourier transformed 
onto the corresponding grid in location space. The r.m.s. strength B is given by 
B 2 = Jq°° dkk 2 (B 2 (k)). We use ub = —11/3, corresponding to Kolmogorov turbulence. 

This turbulent spectrum is applied to three different cases: the first case is an 
all-pervading relatively weak magnetic field ~ InG. For the case of a relatively strong 
magnetic field we performed two simulations: in the first one the field follows a rough 
representation of the Local Supercluster, a pancake profile with scale height of 12 Mpc 
and scale length of 25 Mpc; the center of the profile as at the Virgo cluster, 20 Mpc from 
the observer who is located in the plane of the pancake. The maximal field strength at 
the Virgo cluster equals 0.1/xG. Its strength at the observer is then B = 0.5 x 10~ 7 Gauss. 
In the second simulation we consider an all-pervading field of strength B = 0.1 fiG. 

4. Multi-poles and Autocorrelation function 

For each simulated sky distribution typically 1000 mock data sets consisting of A D b s 
observed events were selected randomly and multiplied with the solid-angle dependent 
exposure function. For each such mock data set or for the real data set we then obtained 
estimators for the spherical harmonic coefficients C (I) and the autocorrelation function 
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Table 2. Experiments considered in the present study. Here, Res. 40 is the angular 
resolution for energies above 40EeV, a^ S ^ is the latitude North (South) of the 
experiment and the angle 9 m is the maximal zenith angle out to which the detector is 
fully efficient. 



Exp Res. 40 a$ a$ 6, 



AGASA 1.6° 
SUGAR 10° 
AUGER 1.° 



+35.5° 



+35° 



45° 

-30.5° 55° 
-39° 60° 



N(6). The estimator for C{1) is defined as 



1 1 j / ^ob s 1 



2 



(5) 



where u>i is the total experimental exposure at arrival direction u l , Af = X)2a s 
is the sum of the weights l/wi, and Yi m {u l ) is the real- valued spherical harmonics 
function taken at direction u l . For a detector at a single site we use the following 
parameterization: 



where a is the latitude of the detector and a m is zero for £ > 1, ir for £ < —1, and 
cos _1 (£) otherwise, where £ = (cos6> m — sin do sin 8)/ [cos do cos 8]. 

In Table El we give a list of the experimental features for the three experiments used 
in the next section. In order to have maximal sky coverage for the large scale multi- 
poles, we included data from the SUGAR array which operated from January 1968 to 
February 1979 in Australia at a latitude of 30.5° South and 149°38' East 34 . 

The estimator for N(9) is defined as 



where S(9) is the solid angle size of the corresponding bin and C = Q e / (A^ b s (iV b s — 
1)), with Q e denoting the solid angle of the sky region where the experiment has non- 
vanishing exposure. In both cases for each / we plot the average over all trials and 
realizations as well as two error bars. The smaller error bar (shown to the left of the 
average) is the statistical error, i.e. the fluctuations due to the finite number N Q b s of 
observed events, averaged over all realizations, while the larger error bar (shown to 
the right of the average) is the "total error" , i.e. the statistical error plus the cosmic 
variance, in other words, the fluctuations due to the finite number of events and the 
variation between different realizations of the magnetic field. 

In the figures shown in Sec.V, the histogram represents the data and the solid line 
represents the analytical prediction for an isotropic distribution in the case of future 



uj(S) oc cos a cos 5 sin a m + a m sin a sin 5 



(6) 
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predictions. Given a set of observed and simulated events we define 



Xn 



E 



as. 



(8) 



where S^data refers to C\ and N(9) defined above, obtained from the real data, and 
S^simu and AS'^simu are the average and standard deviations of these quantities obtained 
from the simulated data sets. This measure of deviation from the average prediction is 
used to obtain an overall likelihood C for the consistency of a given theoretical model 
with an observed data set by counting the fraction of simulated data sets with \n larger 
than the one for the real data. The likelihoods are computed for n=4 in Eq. |HJ 

5. Results 



In the following we compare the results obtained for the simulated UHECR propagation 
scenarios described above with the observational results. As discussed in the previous 
section, the comparison is based on the statistical properties of the simulated 
and observed events, expressed in terms of the angular power spectrum and the 
autocorrelation function of the UHECR arrival distributions. 



0.1000 



0.0100 



0.0010 - 



0.0001 



N obs ( 40. < E/EeV < 1000.)= 99 




Figure 1. The angular power spectrum C(l) as a function of multi-pole I, obtained for 
the AGASA+SUGAR exposure function, for N b s =99 events observed above 40 EcV, 
sampled from 15 simulated configurations with B = 1 nG. The diamonds indicate the 
realization averages, and the left and right error bars represent the statistical and total 
(including cosmic variance due to different realizations) error, see text. The histogram 
represents the AGASA+SUGAR data. The total likelihood for this fit is < 1CT 3 . 



To estimate the true power spectrum from Eq. © requires data with full sky 
coverage and therefore at least two detector sites such as foreseen for the Pierre Auger 
experiment. We therefore combine data from the AG AS A and SUGAR experiments 
which had comparable exposure in the northern and southern hemisphere, respectively. 
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The data set is made of 99 events: the 50 from AGASA (excluding 7 events observed 
by Akeno) and the 49 from Sugar above 4 x 10 19 eV. Usually care has to be taken in 
combining data from two experiments with significantly different angular resolution, see 
Table |21 This is possible only for multi-poles I < 10 which are not sensitive to scales 
< 10°, as pointed out in [HH| . 

In FigQ we compare with the combined AGASA+SUGAR data the angular 
power spectrum predicted by simulations for a scenario with all-pervading fields of 
r.m.s. strength B = 1 nG. The large predicted anisotropy is due to the fact that 
the sources follow the supergalactic plane and that UHECRs do not diffuse on the 
scale of the distance to the dominant, relatively close sources, since their Larmor 
radius r L ~ 110(i?/10 20 el / )(i?/10~ 9 G') _1 Mpc. In addition, the distribution of events 
is concentrated around the directions to a few close sources which give the principal 
contribution to the detected flux. This can be understood by considering, in FigGl the 
number of events detected from source "i" in a given range of energy n^Ei < E < E 2 ), 
normalized to the total number of events detected from all sources in the same energy 
range, n tot (Ei < E < E 2 ), for this weak field case. The labels on the x axis correspond 
to the distances of sources in Mpc whereas at the top of each bin we mark the maximal 
acceleration energy of that source. Almost all particles are detected in the range between 
40-100 EeV, whereas only 12% of the total 75000 particles are detected outside this 
interval. The corresponding histogram also shows that only the sources in the list within 
20 Mpc from Earth give a significant contribution. In case of negligible deflection we 
expect that clusters just reflect the point-like sources; if the number of observed events 
iV b s is larger than the number of contributing sources, as in the present scenario, each 
source contributes on average more than one event and strong clustering is expected. 
In addition, most of the closest sources are within an angular distance less than 5° from 
each other. This implies that most of the events come from the same direction, and 
leads to a predicted autocorrelation function which at scales of a few degrees is about 
30 times the value obtained from AGASA data. 

The contribution of sources injecting above 40 EeV is strongly suppressed below 
this energy by two effects: the closest sources can give a contribution at low energies only 
when their maximal acceleration energy is close to 40 EeV, but they are sparse in the list 
considered here; the source NGC 4342 at 11.4 Mpc with E max = 40 EeV is responsible 
for the first peak in the distribution below 40 EeV in Fig. El The farthest sources in 
contrast do not contribute below ~ 40 EeV because their maximal acceleration energy 
is sufficient to keep the energy of the particles above 40 EeV in the non diffusive regime. 
Fig. 121 shows that the contribution from the farthest sources is actually negligible at all 
energies, due to suppression by a factor proportional to "weight factor /D 2 ", where D is 
the distance to the source. 

So far we have found that the list of sources, as given in Table HJ cannot reproduce 
current data. If typical deflections are small, consistency with the data requires source 
distributions that are more homogeneously distributed than our sample in Table ^ If at 
least some of the sources in our list contribute significantly to the total flux, we would 



Ultra-High Energy Cosmic Rays from Quasar Remnants 



10 



0.50 



W 0.40 



0.30 



0.00 



1 00-1 000 EeV 
10^0 EeV 
40-100 EeV 




Distance(Mpc) 

Figure 2. The number of detected events per source in a given range of energy, 
normalized to the total number of events detected in the same range of energy, averaged 
over 15 realizations, for the simulation with B = 1 nG. The x axis represents the 
distance corresponding to each source in Mpc. At the top of each histogram the 
maximal acceleration energy per source is shown. The solid line corresponds to the 
range (10-40) EeV, the dashed one to the range (40-100) EeV, and the pointed one to 
the range (100-1000) EeV. 



still expect a significant correlation between the arrival direction of the events and the 
positions of some of the objects in our list. 

Following the same approach as in Ref. j2H] we made a simple estimate of the 
positional coincidence between the sources in Table d and the UHECR events observed 
by AGASA with energy above 4 x 10 19 eV. We calculate the number of real coincidences 
in circles centered on the position of a given source. For a circle radius of 4.8°, 
corresponding to a 3a error box in angular resolution, coincidences with the AGASA 
data were only found for two relatively far sources, namely NGC4874 and NGC4889, and 
the source NGC 821 at 24 Mpc, which contribute negligibly to the flux, see Fig. El The 
source with the largest individual contribution (NGC5128) would predict ~ 6 events in 
the range between 40 and 100 EeV, with a Poisson probability for no coincidence of about 
0.1%. This implies that either deflection has to be significant and/or a considerable part 
of the observed flux is due to sources not in our list. We note in this context that hints 
of UHECR arrival direction correlations with another list of quasar remnants have been 
reported in Ref. |2~2~j . 

We now investigate whether stronger magnetic fields, by providing larger angular 
deflection, might provide a better match to the observational data for this same source 
distribution. We performed simulations in turbulent magnetic fields following the profile 
of a sheet centered at 20 Mpc from Earth with strength 0.1/xG at the center of the 
sheet. We recall that infinitely extended profiles are consistent with Faraday rotation 
limits for fields up to fractions of a micro Gauss [24J. In Figs. El and |U we show 
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the autocorrelation function and angular distribution, respectively, above 4 x 10 19 eV, 
predicted by a simulation performed with fields following a pancake profile of scale 
length 25 Mpc and scale height 12 Mpc. The predictions are still not consistent with 
the data for our set of sources. Scenarios for pancake scale heights up to ~ 25 Mpc lead 
to similar results. 



N ob5 ( 40. < E/EeV < 1000.)= 57 




Figure 3. The angular correlation function N(6), eq. Q, as a function of the angular 
distance 6, using a bin size of AO = 1°, obtained for the AGASA exposure function, 
for N D b s =57 events observed above 40 EeV, sampled from 10 simulated realizations for 
a magnetic field following a pancake profile of scale length 25 Mpc and scale height 12 
Mpc with B = Q.lfiG at the center. The diamonds indicate the realization averages, 
and the left and right error bars represent the statistical and total (including cosmic 
variance due to different realizations) errors respectively; see text for explanations. 
The histogram represents the AGASA data. The total likelihood for this fit is 0.01. 



40. < E/EeV < 1000. 




360 180 

right ascension [°] 



Figure 4. Angular distribution in equatorial coordinates on the celestial sphere for 
the same scenario as in Fig. [3] The solid line marks the supergalactic plane. The grey 
scale represents the integral flux per solid angle. The pixel size is 1°; the image has 
been convolved to an angular resolution of 1.6° corresponding to AGASA. 



In the scenario discussed above, the profile of the magnetic field in the sheet is 
such that the field is stronger in the middle plane than on the boundaries of the sheet. 



Ultra-High Energy Cosmic Rays from Quasar Remnants 



12 



In Figs. and El we show the angular distribution and the autocorrelation function, 
respectively, in the case where a magnetic field is distributed homogeneously rather 
than following a profile, although this is not a realistic scenario for our structured 
universe. As a result, when the effects of the boundary are absent, the predictions are 
more consistent with the data and the final distribution appears more isotropic. 

40. < E/EeV < 1000. 




360 180 

right ascension [°] 



Figure 5. Same as FigQ] but for a magnetic field B = 0.1/iG distributed 
homogeneously. 



N obs ( 40. < E/EeV < 1000.)= 57 

1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 r 

15 ■- 




0(deg) 

Figure 6. Same as Fig. but for homogeneous field of 0.1/iG. The total likelihood 
for this fit is 0.18. 

In Fig. [7| we compare the spectrum predicted for this scenario with the AGASA 
and HiRes data. At the highest energies the cut-off is simply due to the fact that the 
maximal injection energy is 257 EeV, whereas we did not attempt to explain the flux 
below ~ 40 EeV as we only included sources accelerating beyond this energy in our 
simulations. Less massive quasar remnants could well significantly contribute in this 
energy range as well. 

The spectrum has been normalized to optimally fit the AGASA data which requires 
an average UHECR emission power of ~ 10 41 erg/s for the nearby sources. At the same 
time the maximum power that can be extracted from a Kerr black hole is given by 
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Figure 7. The realization averaged energy spectra corresponding to Figs. 15161 The 
solid line represents the spectrum that would have been detected by AGASA, using the 
exposure Eq. ©. The dotted lines are the fluctuations due to the different magnetic 
field realizations and the one sigma error bars indicate the AGASA data (full dot) and 
HiRcs data (asterisk). 



L BH ~ 10 40 (a/M) 2 M|S| ergs" 1 [H]. Using Eq. 0, the definition of M and the 
approximation (a/M) ~ 1, as before, we parametrize the UHECR power as 

L CR ~ a Cfl (0.17)10 46 Mf 2 ergs" 1 , (9) 

where acn is the fraction of the power emitted as UHECRs. Comparing these two 
numbers and using typical bulge masses given in Tabled for the closest sources implies 
a required UHECR injection efficiency of about 1% of the accretion rate. 

The spectrum in Fig. shows no very pronounced GZK cut-off. This provides an 
important test for future experiments, especially the Pierre Auger project, which will 
be able to eventually confirm or falsify the presence of a GZK cut-off in the spectrum. 
If this cut-off is not observed, the quasar remnant scenario for UHECR origin here 
investigated could be promising good possibility but with a larger sample of sources 
such as to reproduce the observed isotropy at large scale. Future projects may provide 
sufficient statistics to probe the wiggles seen in Fig. [7| predicted by the monoenergetic 
injection. 

In Fig|Hlwe shown the angular power spectrum for the exposure of the Pierre Auger 
experiment with full sky coverage, assuming 1500 events observed above 40 EeV, for 
the case of an all-pervading 0.1 /iG field. For the exposure function we add Eq. El for 
two sites located at a = —35° and at ao = 39°. The solid line represents the analytical 
prediction for a fully isotropic distribution, predicting Ci ~ (47rA^ r b s )- 1 - This scenario 
predicts an anisotropy that should be easily detectable by the Pierre Auger experiment. 
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Figure 8. The angular power spectrum C(l) as a function of multi-pole I, obtained for 
the two-site Auger exposure function, assuming N o b s =1500 events observed above 40 
EeV, sampled from 10 simulated realizations for an all-pervading turbulent B = 0.1 fj,G 
field. The solid line represents the analytical prediction for an isotropic distribution. 

6. Conclusions 

In the present work we considered a model where the sources of ultra-high energy 
cosmic rays are quasar remnants or dormant AGNs, with underlying supermassive black 
holes as suggested in Ref. [21] ■ We assumed a list of 37 of such objects injecting 
ultra high energy cosmic rays at an energy and with a power determined by their 
mass and accretion properties. We then studied the effects of propagation in different 
extragalactic magnetic fields scenarios on predicted distributions of arrival energies 
and directions. As statistical quantities for this analysis we used spherical multi-poles 
and the autocorrelation function. We found that for a weak magnetic field, of order 
of InG, the predictions appear to be inconsistent with the observed distribution, as 
already pointed out in Ref. [33] , because the magnetic field is too weak to isotropize the 
distribution coming from a limited number of non-uniformly distributed sources, as in 
our list. 

We also found that the contribution from the farthest sources is completely 
negligible even for this weak magnetic field; this is true to an even larger degree for 
stronger fields, consistent with Ref. j33j . 

We found no coincidences between the position of the sources in our list and the 
arrival direction of the AGASA events. In the scenario with a weak magnetic field, 
we thus conclude that the list of quasar remnants considered here can not provide the 
dominant contribution to the observed flux. However, hints of correlations with another 
set of objects subjected to contain quasar remnants have been reported elsewhere [22]. 

These results show that if quasar remnants are the sources of ultra high energy 
cosmic rays, many more of them than considered here must contribute and/or an 
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extended extragalactic magnetic field of ~ 0.1 /iG must exist. We note, however, that 
the number of sources cannot be arbitrarily high, as it would get in conflict with the 
estimates of total density and mass function of supermassive black holes at the current 
epoch pH] . 

We also performed simulations with magnetic fields of ~ 0.1 /zG following pancake 
profiles of various extensions. These two-dimensional sheets were centered at the Virgo 
cluster, at about 20 Mpc from the observer. For a sheet dimension of ~ 20 Mpc in the 
direction towards the observer, and ~ 10 — 20 Mpc perpendicular to it, the predicted 
large scale anisotropy and the autocorrelation beyond a few degrees are too large to be 
consistent with the data. The arrival directions are concentrated along the supergalactic 
plane and the conclusions are basically the same as for the weak field case. Only for all- 
pervading fields of ~ 0.1 yuG with homogeneous properties the predictions become more 
consistent with the data. However, this is unlikely to represent a realistic description 
for our structured universe. 

The quasar remnant scenario tends to predict a relatively hard spectrum with at 
best a mild GZK cut-off which leaves the possibility for this kind of objects to be 
sources at the highest energies. Nevertheless, we recall that these results could be 
affected by the assumption of a mono-energetic injection, which is the more optimistic 
case but not necessarily the most realistic. Furthermore, if the flux at highest energies 
is dominated by just a few sources, mono-energetic or very hard injection spectra may 
lead to conspicuous wiggles in the spectrum. Finally, rough estimates show that about 
1% of the accretion rate emitted in ultra-high energy cosmic rays suffices to reproduce 
the observed flux level. 

As already pointed out in Ref . [2H] , this scenario also has some direct observational 
consequences at much lower energies which could be tested in the near future. The 
emitted spectrum of curvature photons peaks in the TeV band; therefore, these sources 
could be detected by future experiments like GLAST and MAGIC |2B]- At the highest 
energies, the present development of large new detectors, such as the Pierre Auger [To] 
experiment and the space-based air shower detectors such as OWL [Slj and EUSO [38J, 
will considerably increase statistics. All of them are planned to achieve full sky coverage 
and anisotropics predicted by the quasar remnant scenario should be easily detectable. 
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